#include "manager.h"

// calculates: Cij, Dij and nodal volume
// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
// IN:
// OUT:
void calculateCoefficients2D_old(pMesh theMesh, set<int> &setOfDomains, const char *meshFilename)
{
	cout << "calculateCoefficients2D_old\n";
	theMesh->modifyState(1,2);

	// #ifdef MULTIGRID
	theMesh->modifyState(0,1);	// create edge adjacency around nodes
	// #endif

	cout << "Mesh dimension: 2D\n";
	int i,j,numEdges;
	double edgeCenter[3], Dij[2], Cij[3], I[3], J[3], val;

	char *Cij_string[2], *Dij_string[2];
	for (i=0; i<2; i++)
	{
		Cij_string[i] = new char[64];
		Dij_string[i] = new char[64];
	}

    // #ifdef MULTIGRID
    char *deltaDij_string[2];
    for (i=0; i<2; i++)
    {
        deltaDij_string[i] = new char[64];
    }
    // #endif

	pGEntity gEntity;
	vector<pVertex> vertex;
	pEntity edge, face;

	// loop over elements:
	// for each face, calculate all data needed and attach them to the its edges
	// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
	FIter fit = M_faceIter(theMesh);
	while (face = FIter_next(fit))
	{
		gEntity =(pGEntity)F_whatIn((pFace)face);	// get element flag
		setOfDomains.insert(GEN_tag(gEntity));		// store every new domain

		sprintf(Cij_string[0],"Cx-%d",GEN_tag(gEntity));
		sprintf(Cij_string[1],"Cy-%d",GEN_tag(gEntity));

		double faceCenter[3] = {.0, .0, .0};
		F_center(face,faceCenter);								// element centroid

		// loop over all three face's edges to calculate Cij and Dij
		for (i=0; i<3; i++)
		{
			edge = F_edge(face, i);
			edgeCenter[0] = edgeCenter[1] = .0;
			E_center(edge,edgeCenter);				// edge centroid
			M_GetVertices(edge,vertex);
			V_coord(vertex[0],I);
			V_coord(vertex[1],J);

			// edge vector, used as a reference vector
			double IJ[2] = {J[0]-I[0], J[1]-I[1]};

			// vector IJ must point from the smaller vertex ID to the greater
			if ( EN_id(vertex[0]) > EN_id(vertex[1]) )
				for (j=0; j<2; j++) IJ[j] = -IJ[j];

			// vector: from element center to edge middle point, used as a reference vector
			double v[2] = {edgeCenter[0]-faceCenter[0],edgeCenter[1]-faceCenter[1]};

			// Cij is orthogonal to v
			Cij[0] = Cij[1] = .0;
			Cij[0] = v[1];
			Cij[1] = -v[0];

			// Cij must point as if I inside the CV and J outside
			if ( dot(Cij,IJ) <= .0 )
				for (j=0; j<2; j++)	Cij[j] = -Cij[j];

			// associate Cij coefficient to edge
			for (j=0; j<2; j++)
			{
				val = .0;
				EN_getDataDbl(edge,MD_lookupMeshDataId(Cij_string[j]),&val);
				val += Cij[j];
				EN_attachDataDbl(edge,MD_lookupMeshDataId(Cij_string[j]),val);
			}

			// calculate Dij coefficient only for boundary edges.
			// Geometry or domain boundaries
			// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
			pEntity otherFace;
			pGEntity otherFace_gEntity;
			for (j=0; j<E_numFaces(edge); j++)
			{
				otherFace = E_face(edge, j);
				if (face != otherFace){
					otherFace_gEntity =(pGEntity)F_whatIn((pFace)otherFace);
					break;
				}
			}
			// #ifdef MULTIGRID
            // Attaching boundary flag to edge(because of D's Coeff.)
            // and node(because of Priority Queue)
			EN_attachDataInt(edge, MD_lookupMeshDataId("isBoundary"), 0);

            //cout<<"pp-2D: n.Edges: "<<V_numEdges(vertex[0])<<"\n";
            int isBoundary[2] = {0,0};
            // Get the isBoundary from possible previous iteration of the loop
            // through these nodes
            for (int k=0; k<2; k++)
            {
                EN_getDataInt(vertex[k], MD_lookupMeshDataId("isBoundary"), &isBoundary[k]);
            }
            // If it was already defined as boundary(isBoundary = 1), leave it so, otherwise define it as
            // not boundary(isBoundary = 0)
            for (int k=0; k<2; k++)
            {
                if (isBoundary[k] != 1) // Node is not on boundary
                {
                    EN_attachDataInt(vertex[k], MD_lookupMeshDataId("isBoundary"), 0);
                }
            }
            // #endif

			if ( (E_numFaces(edge) == 1) || (GEN_tag(otherFace_gEntity) != GEN_tag(gEntity)) )
			{

                // #ifdef MULTIGRID
				// Attaching boundary flag to edge(because of D's Coeff.)
				// and node(because of Priority Queue)
				EN_attachDataInt(edge, MD_lookupMeshDataId("isBoundary"), 1);
				EN_attachDataInt(vertex[0], MD_lookupMeshDataId("isBoundary"), 1);
				EN_attachDataInt(vertex[1], MD_lookupMeshDataId("isBoundary"), 1);

				// (DEV) - JUST BECAUSE OF AGGLOMERATION, could use a #ifdef...
				// Defining delta D =0, because its median dual, it will be
				// not zero just in the coarse meshes
				sprintf(deltaDij_string[0],"deltaDx-%d",GEN_tag(gEntity));
				sprintf(deltaDij_string[1],"deltaDy-%d",GEN_tag(gEntity));
				for (j=0; j<2; j++)
				{
				    val =.0;
				    EN_attachDataDbl(edge,MD_lookupMeshDataId(deltaDij_string[j]),val);
				}
                // #endif

				Dij[0] = Dij[1] = .0;
				Dij[0] = IJ[1]/2.0;
				Dij[1] = -IJ[0]/2.0;
				// Dij must point to outside element
				if (  dot(Dij,v) <= .0 )
					for (j=0; j<2; j++) Dij[j] = -Dij[j];

				sprintf(Dij_string[0],"Dx-%d",GEN_tag(gEntity));
				sprintf(Dij_string[1],"Dy-%d",GEN_tag(gEntity));

				for (j=0; j<2; j++)
				{
					val = .0;
					EN_getDataDbl(edge,MD_lookupMeshDataId(Dij_string[j]),&val);
					val += Dij[j];
					EN_attachDataDbl(edge,MD_lookupMeshDataId(Dij_string[j]),val);
				}
			}
			vertex.clear();
		}

		// calculate volume of control volume and associate it to elements nodes
		// = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
		char volume_string[64];
		double A = F_area(face)/3.0;							// element area
		M_GetVertices(face,vertex);							// extract face's vertices
		sprintf(volume_string,"volume-%d",GEN_tag(gEntity));
		for (j=0; j<3; j++)
		{
			val = .0;
			EN_getDataDbl(vertex[j],MD_lookupMeshDataId(volume_string),&val);
			val += A;
			EN_attachDataDbl(vertex[j],MD_lookupMeshDataId(volume_string),val);
		}
		vertex.clear();
	}
	FIter_delete(fit);

	// for parallel simulation. Does nothing in serial
	//unifyEdgeCoefficientsMatrices(theMesh);

	for (i=0; i<2; i++)
	{
		delete[] Cij_string[i];
		delete[] Dij_string[i];
	}

	// Export Coefficients
	string str(meshFilename);
	string::size_type start = str.find_last_of("/");
	start = (start == string::npos)?0:start;
	string::size_type end = str.find_last_of(".");
	char buffer[256];
	memset( buffer, '\0', 256 );
	str.copy(buffer,end-start,start);
	char filename[256];
	sprintf(filename,"preprocessor-datafiles/%s.dat",buffer);
	exportCoefficients(theMesh,2,setOfDomains,filename);

	// Calling the agglomeration procedure
	// (DEV) - We must have some conditional defining if the user
	// wishes to use the agglomeration procedure or not
	outputPrimal(theMesh, 2, filename,-1, theMesh);
	agg2D(theMesh, setOfDomains, meshFilename);
}
